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ABSTRACT 

We present a hybrid technique of N-body simulation to deal with coUisionless stellar 
systems having an inhomogeneous global structure. We combine a treecode and a self- 
consistent field code such that each of the codes model a different component of the 
system being investigated. The treecode is suited to treatment of dynamically cold 
or clumpy systems which may undergo significant evolution within a dynamically hot 
system. The hot system is appropriately evolved by the self-consistent field code. This 
combined code is particularly suited to a number of problems in galactic dynamics. 
Applications of the code to these problems are briefly discussed. 

Key words: methods: numerical - stellar dynamics ~ galaxies: kinematics and dy- 
namics - galaxies: interactions 



1 INTRODUCTION 

We concern ourselves with the general dynamical evolution 
of a se lf-gravitatinK stellar system dominated to the first- 



order bv a stable or slowly evolving spheriodal mass distribu- 



tion but with significant and dynamically distinct substruc- 
ture. We can refer to this as an inhomogeneous spheroidal 
system, e.g. a disk or other structure embedded within a 
dark matter halo, an encounter between an elliptical galaxy 
and less-massive companion, sinking satellites, rings, and 
fine structure in ellipticals. The aim in this work is to de- 
velop an efficient and handy method of modelling the com- 
plex dynamical evolution of the types of systems mentioned, 
given the resources which are commonly at hand. 

An established method of investigating the dynamical 
evolution of stellar systems is through the use of N-body sim- 
ulations. Such simulations have become increasingly sophis- 
ticated over the last decade, enabling detailed experiments 
to be performed on N-body models of stellar systems. Im- 
provements have been a result of greater computing power 
and more cunning algorithms for following the time evolu- 
tion of a system of particles. A major consideration for the 
researcher is to match the computing resources available to 
an N-body method suited to modelling the system under 
consideration. 

The majority of numerical treatments of stellar dynam- 
ics rely upon the assumption that stellar systems on the 
scale of galaxies are collision free. This is based on the fact 
that the two-body relaxation time of a star, ireiax, is many 
magnitudes larger than the age of the galaxy which contains 
it. The appropriate globally averaged estimate is given by. 
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for a system of TV bodies where the crossing tim e, tc 



the time for a particle to cross the system once (Binney & 



Tremaine 1987) 



The two-body relaxation rate of a body in self- 
gravitating systems depends in part on the local density 
and velocity dispersion. Relaxation of the orbits of individ- 
ual particles will occur more quickly in regions of higher den- 
sity or lower velocity dispersion. Other authors present more 
detailed discussions of relaxation processes {e.g. Farouki & 
Salpeter 1994; Huang et al. 1993). If one is concentrating on 
short timescale evolution in a region of fine substructure in 
an otherwise dynamically stable or only slowly evolving sys- 
tem, one may find that the resolution locally is insufficient 
to accurately describe the detailed evolution. This provides 
the motivation to combine simulation techniques which will 
deal seperately but efficiently with different components of 
an N-body system. 



1.1 The N-body Problem 

The dynamical evolution of a system of coUisionless particles 
is described by the coUisionless Boltzmann equation. 



and Poisson's equation. 



(2) 



(3) 



where / is the distribution function of the particles, describ- 
ing their position r and velocity t; at a time i, and $ is the 
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gravitational potential at (r,t) due to all the particles. The 
density of the system is related to / by 



fd^v. 



mizing statistical fluctuations in particle distribution. Ifow- 
ever spatial resolution is constrained by the grid spacing. 
A hybrid technique has been developed incorporating PP 



M. 



and PM schemes, unsur prisingly referred to as P M (East 



wood & Hockney 1974). This technique has been usefully 



G is the universal gravitational constant and hereafter is 
taken t o be unity. The acceleration can then be found from 



the potsntial 



a = -V$. 



(5) 



Apart from a number of exact equilibrium solutions of 
the coUisionless Boltzmann equations, in general this func- 
tion of seven independent variables cannot be solved. A fea- 
sible way of solving for the time evolution of equations (pi) 
and (m) is to construct an N-body realisation of the system 
by sampling the phase space {r,v) N times, subject to the 
probability density /(r, v). The N-body system of particles 
is then evolved according to Newton's laws. It is desirable 
for the sample A^ to be as large as possible to reduce the ef- 
fects of statistical noise in the sample, lessening the effects of 
numerical two-body relaxation, and increasing the possible 
spatial resolution. Memory and processing time of comput- 
ing res ources constrain the value of N that is achievable. 



and th us N-body codes generally employ various approxi- 
mations to counter the problems raised with smaller N . In 
general techniques address either one problem or another, 
and so choice of appropriate algorithms is very important. 
One must ensure that the limitations of any particular al- 
gorithm does not invalidate it's application to the physical 
system under consideration. 

In the following section we outline some of the major 
techniques that have been employed in the study of stellar 
dynamical problems [see also Hernquist (1987) for a com- 
prehensive review], and explain why it is we implemented 
the two methods used in SCFTREE. 



1.2 Techniques 

Perhaps the most straightforward way to evolve the system 
is to calculate directly all the interparticle accelerations. The 
combined acceleration, a^, on a particle is 



JV 

E 
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[|r, -r,|2 + e2]3/2' 



(6) 



where r^, rj, rrii, and rrij are the positions and masses of 
the particles i and j, e is the softening parameter. These 
particle-particle (PP) or direct summation methods have 
two important aspects in their favour, the resolution scale 
is determined solely by e, and there are no constraints on 
the geometry of the system. The greatest drawback is in 
computational cost, at best the CPU time scales as 0(N ^ 



Integ r ation techniqu es have become quite efficient (Aarseth 
1973 lAarseth 19851), but CPU expensive for A^ ^ 10* (al- 



though new dedicated hardware has improved the situation, 
see below). 

An extensio n of the PP technique is the particle-mesh 
(PM) technique ([Hockney fc Eastwood 198l|). This method 



imposes a grid upon the system and densities are a ssigned 



to each grid cell . Applying Fast Fourier Transform (Cooley 



&i Turkey 1965) to PM methods makes them highly effi- 
cient at dealing with large numbers of particles, thus mini- 



implemented for simulations of large-scale structure, where 



high density contrasts are expected (Efstathiou & East 



wood 1981), however when densities get high the number 



of close neighbours, which are dealt with by the PP algo- 
rithm, becomes large and prohibitively lengthen the compu- 
tation time. Geometric constraints are also imposed by the 
existence of the fixed grid. Further refinement to the P^M 
technique has been to introduce an adaptive grid (AP'^M, 
Couchman 1991, Couchman et al. 1995). Examples of fur- 
ther variants of the adaptive mesh technique are the adap- 
tive Particle-Muhiple-Mesh (PM^) of Gelato, Chernoff, & 
Wasserman (1996), and the hydrodynamic and N-body un- 
structured adaptive mesh of Xu (1995). 

For systems with a fairly high degree of symmetry it is 
possible to represent the potential of the system as a series of 
terms in a multipole expansion about the centre of symme- 
try of the system. V arious basis functions for the expansion 
have been employed ( Clutton-Brock 197S ;_ van Albada fc van 



Gorkom 1977; McGlynn 1984; Hernquist fc Ostriker 1992) 



depending on the global geometry of the system being mod- 
elled. The expansion is truncated at a specified order, n, 
this governs the effective resolution of the technique. The 
primary advantage of this technique is that computation 
time scales as O( nN), making large A'' an attractive pos- 
sibility. Weinberg (1996) has recently presented an advance 



on the expansion technique whereby the expansion basis and 
number of expansion terms are matched to the system dur- 
ing its time-evolution. The main drawback is that expansion 
techniques do not consider individual particle-particle inter- 
actions and so local substructure is ostensively suppressed. 
The SCF method of Hernquist & Ostriker (1992) (hereafter 
HO) , is particularly suited to the observed mass distribution 
in ellipt ical s, we expand upon this method in some more de- 
tail in Q. 

The advantages of the particle in a field approach of 
expansion techniques and the geometrical fiexibility of PP 
methods are combined in what are generically described as 
tree codes. At close range PP interaction forces are cal- 
culated explicitly, but as the range increases, particles are 
grouped together in larger and larger clusters and their com- 
bined potentials are approximated by truncated expansions. 
For each individual particle the force on a particle is the 
sum of progressively more distant particle-cluster interac- 
tions. These methods are known as tree methods due to the 
division of the particle data into successively smaller and 
smaller clusters until a single particle is reached. The ad- 
vantage behind the method is that the time scales involved 
scale as C'(A''log A''), providing significant improvements in 
efficiency over PP methods. There have been several meth- 
ods of constructin g the tree structure from particle data [e.g . 
the AJP method ( Appel 1981] ; pfernigan 198^; [Porter 1985 | ) , 
the BH hierarch ical t hree method ( Barnes & H ut 198 6 ), see 
also Hernquist ([l987[), and Warren fc Salmon ([199*^. The 



BH method has proved to be the most popular, its method 
of tree construction is more organised and efficient, although 
not necessarily as accurate as the AJ P a pproach. The hier- 



archical tree method is described in 52.1 
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Recent modifications (Warren & Salmon 1993) to the 



hierarchical tree process has given greater control over the 
accuracy of the force calculations for the same computa- 
tional cost of C»(iVlogiV). See also Salmon & Warren ( |l994| ) 
for detailed discussion on error analysis of treecodes. 

A potentially promisin g approach has been the Fast 
Multipole Method (FMM) (iGreengard & Rokhlin 19871). In 



principle the FMM is similar to other tree methods described 
above, but in calculating the potential at a point, the algo- 
rithm passes from the root node to the leaves, accumulating 
information from the cells at the same level according to a 
range criterion. This method has been successful as used in 
two-dimensional, unclustered applications, but h as failed to 
be as efficient in three-dimensional applications ( Schmidt & 
Lee 19d^L 



Re :ent advances in computer architecture have seen a 
shift to massively parallel machines. Simply put, the com- 
puting load is split between A'" separate processors simul- 
taneously, increasing the computational speed by a factor 
of A'^. The method by which it is split varies depending on 
the algorithms implemented, some being more suitable tha n 
others, such as the SCF method (Hernquist et al. 1995). 
Although intrinsically more tricky to parallelise, tree codes 
have al so been successfully implemented on parallel ma- 
chines ([Salmon 199l|; |01sen fc Dorband 1994 [Warren 1994 
Dubinski 199^ ; pavi 19971 ). 

An exciting prospect in N-body studies has been the 
introduction of special purpose computer chips which are 
dedicated to the force calculations. In the simplest case, for 
example in a PP code, one would supplement calls to the ac- 
celeration subroutine directly by calls to the chip. Details of 
the range of "GRAPE" and "HARP" boar ds tha t have been 
developed can be found in Ebisuzaki et al. (1993). These de- 



velopments have also been implemented in codes based on 
Tree algorithms {e.g. Steinmetz 1996). 

For the generic problem of a spheroidal system with 
fine substructure, the various advantages and disadvantages 
exhibited by each of the methods point towards a combi- 
nation of the SCF method and BH treecode scheme. These 
have the advantage of having been widely tested in a num- 
ber of astrophysical applications and have been successfully 
implemented on commonly available unix workstations. In 
§e| we outline in more detail the principles behind the spe- 
cific Tree and SCF methods used here. Following this we 
describe how these two separate algorithms were combined 
to form the new hybrid technique, and how we implemented 
this new SCFTREE code. We thoroughly test the efficiency 
and accuracy of SCFTREE, and present the results in §pl Fi- 
nally in §y we discuss the range of applicability of this new 
technique to a wide variety of dynamical problems. 



2 PUTTING TREES IN FIELDS 



2.1 The Hierarchical Tree Code 



The nature of the hierarchical tree method (Barnes & Hut 



1986 ), as introduced in the previous section lends itself to 
efficient and logical coding of its algorithm. 

The basic principle in improving the direct summation 
(PP) technique is to group together long range interactions. 
In the BH tree code particles are contained within a cubic 



volume, known as the root cell. The tree code algorithm gets 
its name from the way the particles are grouped together in a 
hierarchical level of cells, with the root cell being at the top. 
This root cell is subdivided into eight further cells, the next 
level down in the hierarchy of cells. If any of these cells con- 
tain more than one particle, that cell is further subdivided 
into eight. This 'tree building' process continues until the 
subdivision of cells at a particular level in the tree building 
can go no further, i.e. the final hierarchical tree level is com- 
posed of cells which contain only one particle. For all cells 
on every level which contain more than one particle, a mul- 
tipole expansion is performed about its centre of mass. This 
expansion is typically truncated at the quadrupole term, al- 
though some implementations of the t ree algorithm include 
highe r terms, e.g. the octupole term (McMillan & Aarseth 



199J). Here we truncate the expansion at the quadrupole 
term, so that for each cell we calculate its centre of mass 



and its quadrupole moment, Q, (Goldstein 1980). The ac- 
celeration at a position Tc from the centre of mass of the 
cell can then be expressed as 



ai — G 



McTc , Q ■ Tc 



5(r, 
2 



Q r. 



(7) 



For any particle in the system it, in effect, now sees a hier- 
archy of cells of different sizes and different distances. If s is 
the size of a cell and d is the distance from the particle to 
the cell, the simple criterion, 



> 



d' 



(8) 



where 6 is known as the tolerance parameter, can be ap- 
plied to decide whether the interaction between particles in 
that cell and the one particle is adequately approximated by 
the acceleration given by equationdTI), or whether the par- 
ticle should 'look' at the lower level of smaller cells. If the 
cell contains a single particle then the acceleration is simply 
given by equation (0). The tolerance parameter, 6, can be 
set by the investigator to control the accuracy of the approx- 
imation. The usual process of calculating forces on a particle 
is to 'walk' through the tree from the root cell downwards 
forming a list of valid interactions between cells and the par- 
ticle. The length of this list is proportional to log N , hence 
the CPU time for a single timestep scales with 0{N\o^N). 
This is what makes the tree method a distinct improvement 
over the PP method. 



2.2 Self-Consistent Field Code 

As far as a star in a galaxy is concerned, the dominant com- 
ponent of the gravitational potential which it experiences 
comes from the averaged field of billions of distant gravita- 
tional sources. It may appear rather odd that on the whole, 
the simulation of the gravitational dynamics of galaxies em- 
ploys a method of basically summing up all the individ- 
ual particle-particle interactions, each one of which is j ust 
about negligible. The expansion techniques described in §1.2 
appear more natural, as the force calculation on each par- 
ticle considers the effect of the mean gravitational field of 
the whole system directly acting on the particle. One might 
expect that the effects of two-body relaxation would be os- 
tensively mitigated. On the other hand the model system 
still consists of many fewer bodies than a real system, and 
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SCF PARTICLE 




SCF system 



TREE PARTICLE 



Figure 1. Schematic representation of tlie interaction between the SCF system and Tree system. The potentials $i and <1>4 are d ue to 
the expanded SCF potential, see equation (12), at a position relative to the centre of mass of the SCF system (as defi ned in §2.3). The 
potentials $2 and <J>3 are due to the tree system of particles. The calculation of the tree potential is described in §p.]|. 



thus the statistical fluctuations brought about by the finite 
number of particles will cause the expanded gravitational 
field to fluctuate and thus the overall effect on the system 
approaches that of one that has experienced a correspond- 
ing amount of two-body relaxation. The major advantage of 
the SCF method is that the speed of the computation scales 
as 0{N). Therefore it is possible to increase the practical 
limit on the number of particles, A'', and thus any statistical 
fluctuations in the expanded potential are reduced. 

The principle behind the SCF approach is to represent 
the potential of the system as a truncated series of terms 
of an expanded set of basis functions. It is possible to find 
the coefficients of such an expansion from the known den- 
sity field 'sampled' by the particles. Poisson's equation is 
solved for the set of basis functions, and coefficients of the 
expanded potential can be found. From the potential ex- 
pansion the acceleration of any particle can be directly cal- 
culated. This technique had been applied in vario us guises 



and va rious d egrees to a limited extent in the pas t (Glutton 



Brock [973; van Albada & van Gorkom 1977; McGlynn 



1984). "?he philosophy behind the SCF method of Hernquist 



& Ostriker (1992) was to use a set of basis function where 
the lowest order terms well represent the observed distribu- 
tion in spheroidal galaxies. This method has subsequently 
been successfully utilised in further studies, [e.g. Joh nston , 
Spergel, & Hernquist (1995); Hoz umi fc Hernquist ( 1995 ), 
Hernquist, Sigurds son, a nd Bryan ( 1995|) , Sigurdsson, Hern- 
quist, & Quinlan (1199E|)1. Below we present a quick resume 
of the method of ([Hernquist fc Ostriker 1992 ) . 



Hernquist (1990) demonstrated that a density-potential 
pair exists such that their properties provide a very good 



approximation to the actual properties of observed galaxies, 
i.e. the R^'"^ law. They are: 

Ma 1 



p(r) 



and 



(r) 



2n r (r + a)^ ' 



M 



r + a 



(9) 



(10) 



The total mass is M, and a is the scale length related to 
the half-mass radius such that ri/2 = (1 -I- \/2)a. A great 
advantage of this density-potential pair is that many of its 
properties can be derived analytically. 

A biorthogonal basis set is constructed for the density 
and the potential, and they can be written as expanded se- 



Pir) = y^^AnlmPnlmjr), 



and 



(11) 



(12) 



where n, I, and m are equivalent to "quantum" numbers in 
the expansion, one radial, and two aiigular. Each pair pnim 
and ^nim satisfy Possion's equation (|3|). 

The lowest order terms are set to be those of the as- 



sumed model, i.e. 
where G — M = 



Pooo 
a 



and $000 = -TTz:, 



2tt r (1 + r)^ """ " """ ~ 1 + r ' 

= 1. It is possible to construct the 
higher order terms of the expanded potential-density series 
'^■nim.ij^ s-nd p„im{r), and to derive an expression for cal- 
culation of the coefficients A„im ■ The interested reader may 
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r efer t o the clear presentation given in Hernquist & Ostriker 
( [L992| ). 

Once all the Anim have been calculated from the known 
coordinates of all the particles, the potential at the loca- 
tion of any particle can be evaluated using equation (|l2|), 
and thus the acceleration of the particle can be found using 
equation (B. 

For the purposes of our hybrid code, the basis set de- 
rived by Hernquist & Ostriker (1992) is an ideal choice. 
Other basis sets may be more useful for other systems such 
as discs, or objects that are not well-approximated by the 
7?^'* profile. [Other basis sets are discussed elsewhere e.g. 
Saha (1993); Earn (1995); Earn & Sellwood (1995); Zhao 
(1996); Syer (1995)] 



2.3 Combining Tree and SCF 

The problem we are faced with is how to combine the two 
codes. It is a relatively simple task, all we need to achieve is 



''(on tree partiele) 



''(due to Tree) 



''(due to SCF expansion) 5 



(13) 



fl(on SCF partiele) — Ifduc to SCF expansion) + fl(due to Tree) • (14) 

The left-hand sides of equations {UM and (|l4) are the total 
accelerations on a particle in the TREE and SCF systems 
respectively. The first terms on the right hand side are the 
accelerations on the particles due to their own systems, the 
terms on the far right are the extra accelerations due to 
the other system of particles. This is schematically shown in 
Figure (g). 

When implementing this there are a couple of points to 
bear in mind. Firstly we need to take account of the response 
of the SCF system as a whole to the TREE system, i.e. 
the centre of mass will not be fixed. So we have to make 
sure that the SCF expansion is calculated about the centre 
of mass, or symmetry, of the SCF system of particles. As 
the SCF system evolves it may become asymmetric to some 
degree, hence we need to ensure that the expansion takes 
place about the most tightly bound region. We achieve this 
by calculating the centre of mass only from the contribution 
of particles which have less than the average particle energy. 
This is a parameter which may affect the accuracy of the 
expansion. It is possible to vary the energy level cut o ff for 
calculation of centre of mass. We discuss this further in 53.5 



Our concern now is how to treat the Tree system of par- 
ticles with respect to the individual SCF particles. There 
are two options. Firstly, we may perform an expansion of 
the tree system much in the same way as the SCF system, 
then each SCF particle will view the potential due to the 
Tree system as a truncated set of terms in an expanded po- 
tential. For encounters that are non-penetrating this may 
prove a quick and reasonable approximation, but because 
the tree system is a priori not assumed to have any particu- 
lar geometry, and likely to have a widely disrupted one, we 
would not expect it to be accurately represented in just a 
few terms of an expanded potential. Since we also desire our 
models to describe encounters and interacting systems we 
have to be aware that particle-particle interactions would 
not be accounted for if we used an SCF expanded Tree po- 
tential, e.g. the effects of dynamical friction would be poorly 
represented. Our other, more instinctive, option is to treat 
the system of Tree particles simply as a tree from the point 



of view of the SCF particles. This is a much more logical 
procedure as the tree has already been constructed and all 
that needs to be done for each SCF particle is to construct 
a tree interaction list and then sum the over particle-cell 
forces in the list. 

The main calculations performed by the code in each 
time step of the system are outlined in the following sections. 



2. 3. 1 Expansion Centre of the SCF System 

As stated above the centre of mass of the SCF system is 
updated every timestep. Its position is calculated from the 
most tightly bound particles and is subsequently used as the 
origin of the SCF expansion when producing the coefficients. 



2.3.2 Tree Particle - SCF Interaction 

The position of each tree particle, r -^ (r, 6, 0), with respect 
to the SCF expansion centre, is passed to a routine which 
calculates $(r), i.e. equation (|l3). The Anim have been al- 
ready found for the case of the SCF system. The resulting 
acceleration of the Tree particle due to the SCF system is 
subsequently calculated. 



2.3.3 Acceleration of Tree Particles due to Tree 



As described in §2.1, a tree data structure is built from the 
particles in the tree system. Then for each particle an list of 
interactions between itself and the tree cells is formed sub- 
ject to the tolerance parameter, 6. Finally the cell-particle 
accelerations for each element in the list are calculated, 
summed, and added to the SCF acceleration from 5 2.3.2. 



2.3.i Acceleration of SCF Particles due to SCF 

As described in §p^ the coefficients, Anim, of the SCF ex- 
pansion are calculated from the positions of the SCF parti- 
cles. The acceleration due to the SCF system on each SCF 
particle is then found by applications of equations (121) and 



2.3.5 SCF Particle - Tree Interaction 

For each SCF particle we build a tree interaction list. Each 
level of tree cells from the largest (root) to the smallest (sin- 
gle particle) is examined. If the criterion of inequality (H) is 
satisfied then the cell is added to the interaction list, other- 
wise the next level down is examined. Once the interaction 
list is formed, the list is looped through, summing the ac- 
celerations between the SCF particle and the tree cell as 
calculated by equation (m , and adding this to the accelera- 
tion due to the SCF particles themselves. 



2. 3. 6 summary 

For each timestep: 

(1) Find centre of mass of SCF system to use as centre 
of expansion. 



© 1996 RAS, MNRAS 000, 0-|I| 



6 Vine and Sigurdsson 



(2) Calculate the coefficients Anim of tlie terms in tlie 

SCF expansion. 

(3) Build the Tree from all the Tree particles 

(4) For each Tree particle: 

a) Calculate '^scpirTr^e)- 

b) Form Tree interaction list between TREE particle 

and Tree system subject to 6. 

C) Calculate '^Tree(l'Tree)- 

d) Calculate acceleration of Tree particle. 

(5) For each SCF particle: 

a) Calculate ^scFirscp)- 

b) Form Tree interaction list between SCF particle 

and Tree subject to 6. 

C) Calculate $Tree(r'SCF)- 

d) Calculate acceleration of SCF particle. 

(6) Update positions and velocities of all particles. 

2.4 Implementation 

The aim in this project was to combine two separate codes 
into one hybrid code. The process was made much easier 
by the availability to us of optimised versions of both the 
treecode and SCF code written in Fortran 77. We acknowl- 
edge the author, Lars Hernquist, for making these codes 
freely available. 

The two codes were made compatible and carefully 
purged of any overlapping routines and variables ensuring 
no cause for confusion existed. The extra routines added 
were: calculation of the SCF system centre of mass correc- 
tion; formation of tree interaction list for an SCF particle; 
calculation of the acceleration on a SCF particle from the 
tree list; calculation of the acceleration on a Tree particle 
from the SCF expanded potential; and a number of modifi- 
cations to the output information. 

The final version of the SCFTREE code was imple- 
mented on a Sun Sparc20 workstation running Solaris 2.4, 
and being written in standard F77 should be fully portable. 



3 TESTS 

Having two sets of particles that are evolved under two dif- 
ferent numerical schemes one has to be very careful that the 
system as a whole is behaving as a realistic model of the 
system it represents. We conduct a variety of stringent tests 
on accuracy and efficiency which display the applicability 
of this code to the systems it weis designed to deal with. In 
this section we quantify the efficiency of the code, check its 
validity with respect to the constants of motion, and put its 
performance to the test in a number of dynamical examples. 



3.1 Models 

Three distributions of particles are realised in order to test 
the stability of spherical systems and compare the ability 
of SCFTREE to handle distributions departing from the 
Hernquist distribution (the SCF basi s set). The d istribu- 



halos (Kuijken & Dubinski 1994), and as such is our pre- 
ferred model in many of the tests. With reference to the 
parameters in Kuijken & Dubinski (1994) we use the follow- 
ing for all our Lowered Evans models: 'I'o = —5.0, vo — 1.5, 
g = 1.0, {r,/rkf = 1.0, ra = 1.0. 

In these tests the models are populated with up to 10^ 
particles and a specified fraction of them are designated as 
particles to be dealt with by the Tree part of the code, re- 
ferred to hereafter as 'Tree particles'. The remainder of the 
particles are dealt with by the SCF part of the code and 
referred to as 'SCF particles'. 

A nalog ous to the investigations of Hozumi and Hern- 
quist (1995) who perform tests on the pure SCF code, we 
examine the behaviour of SCFTREE in systems which are 
not in equilibrium. For this purpose we construct a uniform 
density sphere with velocity dispersions assigned according 
to a specified initial virial ratio, the evolution of the system 
is followed until it reaches its final relaxed state. 

Extending these tests of SCFTREE to more interesting 
systems we perform test runs on a disk system. A disk, bulge, 
and halo model is set up as prescribed in Dubinski and Kui- 
jken (1995). Such a system is weU suited to the SCFTREE 
technique whereby the disk particles are assigned to the Tree 
and the halo and bulge particles are assigned to the SCF 
system. Here we merely examine the effect of increasing the 
number of halo particles on the stability of the disk. 



3.2 Timing 

The individual timing performances of tree and SCF codes 
have been shown to scale t o O(Ntree^oRNtree), (Barnes & 
Hut 1986| ), and 0{Nscf), ( [Hernquist fc Ostriker 1992| ), re- 
spectively. For the combined code we have to consider the 
CPU time spent on the interaction between the two compo- 
nents. The calculation of the Tree force on each SCF particle 
would, at worst, scale as 0{Nscf log Ntree), and that of the 
SCF force on the tree particles scales as 0{Ntree). So overall 
we would expect the CPU time to scale as 



O [NsCf{1 + logNtree) + Ntre^C^ + log NsCf)] ■ 



(15) 



tions are the Hernquist de nsity profile ( Hernquist 199C ) , the 
Plummer dens ity profile ([Plummer 1911 ), and the Lowered 



Evans model ( Evans 1993[ ). The Lowered Evans model has 
been shown to be useful for the practical modelling of dark 



As the major advantage of SCFTREE is expected to be its 
efficiency in dealing with systems with substructure, we per- 
form tests varying the fraction of tree particles in the system 
and changing the initial distribution of these particles. The 
models used here for the performance checks are the Low- 
ered Evans models. 

We examine how the timing of SCFTREE behaves with 
the total number of particles. The fraction of tree particles 
is set at 10% of the total and they are distributed randomly 
throughout the system. The random distribution of parti- 
cles implies that there will not be any advantage in time 
saved due to the Tree structure. This is demonstrated in 
Figure 2(a) where we compare purely Tree, purely SCF, 
and SCFTREE timings as the total number of particles is 
increased up to 10^. As expected no advantage is offered 
by SCFTREE with randomly distributed particles over the 
pure Tree treatment. 

The great advantage of SCFTREE is demonstrated in 
Figure 2(b) where the Tree particles are assigned to a satel- 
lite system orbiting a more massive system composed of the 
SCF particles. This is compared to a pure Tree treatment 
of the same binary system. The Figure clearly demonstrates 
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Figure 2. CPU per timestep performances for SCFTREE. 

(a) SCFTREE is compared to pure Tree and pure SCF codes, in the SCFTREE runs 10% of the particles are designated as Tree particles 
and distributed randomly through the system (A lowered Evans model with 20000 particles, dt=0.05). 

(b) Binary system of spherical galaxies in a circular orbit. Mass ratio 10:1. In the SCFTREE runs the smaller satellite is composed of 
only Tree particles and the larger body composed of only SCF Particles. 

(c) Variation of the fraction of Tree particles, Ntrcc/NscF distributed randomly in a Lowered Evans model. 




Figure 3. The magnitude of 2-body relaxation in the system 
versus time. The relaxation measure we use is (\AE/E\'^\, defined 
in equation (hq) . The model is a Lowered Evans model containing 
20,000 particles, 10% being randomly distributed tree particles. 




time 

Figure 4. An example of the total fractional energy change, {Et — 
Eo)/ Eq, as a function of time, for a Lowered Evans model, 20000 
particles with 10& in Tree. SCF truncation parameters are ra = 16 
and / = 6. SCFTREE conserves energy to within duration of this 



the strength of SCFTREE when the Tree particles are as- 
signed to a distinct substructure. We see that for SCFTREE 
the CPU time per timestep scales approximately as 0{N), 
a substantial improvement over the pure Tree timings. 

Finally we examine the performance of the code as the 
fraction of Tree particles is increased. Figure 2(c) shows the 
CPU time spent both on the whole SCFTREE timestep and 
just on the Tree timestep in a system with a randomly dis- 
tributed fraction of Tree particles. It can be seen that for 
such a system the fraction of Tree particles must remain 
below ~ 10% to gain any reasonable advantages. 



3.3 Relcixation 

It is instructive to follow the behaviour of two-body relax- 
ation as the system evolves. As an indicator of two-body re- 
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laxation we take the mean square fractional energy change 
of the particles, R2-body = (\AE/E\''^\. We define this as 



R2-body 



N 

-y 



-Ei,t — -El 



i,0 



Ei. 



(16) 



Ei,t and Ei^ are the total energy (kinetic plus potential) 
of particle i at times t and respectively. In Figure (|3|) 
we see that by our measure the average relaxation remains 
less than 1% for the period of the simulation. The initial 
sharp rise is expected due to the initial random placing of 
the particles, meaning that many pairs of particles will have 
initial separations less than the tree softening length, e. 
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Figure 5. How the individual particle energies behave with respect to different initial equilibrium models. Results are calculated 
t = 10 to t = 50. All models consist of a total of 20000 particles, 2000 of which arc randomly assigned to the Tree system. The 
expansion is truncated at n = 6 and 1 = 2. 
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Figure 6. Showing the effect on individual particle errors as the fraction of Tree particles in the simulation is reduced. This demonstrates 
that it is the presence of the Tree particles in the system which cause the vertical "bands" associated with certain energies. The relative 
energy difference per particle, dE, is calculated from t = 10 to t = 50. These models are Lowered Evans models with a total of 20000 
particles. 



3.4 Energy Conservation 

Relaxation tells us how close our models approach the be- 
haviour of a collisionless system, providing us with a quan- 
titive estimate of over what timescales our simulations can 
be applied. 

We may also look at the change in the energy of each 



particle over a suitable time period and compare the be- 
haviour of SCF and tree particles. This gives another gauge 
of relaxation of particles this time as a function of the bind- 
ing energy, E. We follow a similar analysis as H092 and plot 
AE/E vs E for both SCF and tree particles. Figured^) shows 
the results for the 3 spheroidal models: Hernquist; Plummer; 
and Lowered Evans. Each model is populated with 20000 
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particles and the energy change is fractional difference of 
the particle energy at times t and to- One will notice all 
three SCF particle plots exhibit a vertical band structure at 
certain binding energies. The cause of some particles being 
more liable to undergo two-body relaxation than other is 
made clearer by Figure(H) which shows the same plots for 
three Lowered Evans models with 0, 500, and 1000 tree par- 
ticles respectively. The greater the fraction of tree particles 
the wider and more distinct the vertical features that appear 
in the regions of low binding energies. The first point to bear 
in mind is that at low energies the fractional energy changes 
will become large because the binding energy is approaching 
zero. Secondly it is clear that greater the number of tree par- 
ticles the greater the fractional energy change. This is due to 
the fact that individual SCF particles interact directly with 
the tree particles, and so the more tree particles the greater 
the energy exchange between the two components. Finally 
the bands that are seen appear as a result of truncating 
the SCF expansion. The expanded potential of the system 
is derived from the actual positions of all the constituent 
SCF particles. The actual binding energy of a particle may 
not agree exactly with the potential given by the expansion, 
causing an error in the resulting acceleration. There will be 
certain values of binding energy which correspond to the di- 
vergences in the truncated expanded potential from the true 
value, thus the particles at those binding energies will expe- 
rience greater errors in acceleration. Clearly the presence of 
Tree particles in the system exacerbates this effect. However 
because these effects are symmetric the overall global errors 
remain minimal. 

This is confirmed by monitoring the total energy of the 
system. FigureM) traces the fractional variation in total en- 
ergy of the system, {Et — EtQ)/EtQ. The model is a Lowered 
Evans with 20000 particles, 10% being Tree particles, and 
timestep dt — 0.05. Energy is conserved in this system to 
within 0.8% over a period of 50 time units. 



3.5 Momenta 

As with all expansion codes angular and linear momentum 
are intrinsically not conserved exactly, due to approxima- 
tions in the force calculation. We might expect it to be more 
pronounced in the SCFTREE case where the particles in the 
Tree and SCF codes do not respond to each other equally 
and oppositely. Linear momentum is known not to be con- 
served in the pure SCF case, and can be taken care of by 
recentering the particles. In this case, every expansion of 
the SCF system is taken about its centre of mass, alleviat- 
ing the need to recentre the particle. We check the accuracy 
of this modification by monitoring the conservation of linear 
momentum of the SCF centre of mass of a spherical stellar 
system, and the separation of a binary system of spherical 
galaxies (one SCF; one Tree) in a circular and stable or- 
bit. Figure nO shows the separation between an SCF and a 
Tree spheroid of equal mass, both with 2500 particles, in a 
circular orbit. 

We show an example of the evolution of the total angu- 
lar momentum of the system in Figure (0a) and although |L 
is intrinsically near zero we see that it varies no more than 
1% over the length of the run. Figure (iTp) shows the abso- 
lute variation of the ^-component of L. The system in this 




Figure 7. Plot (a) shows the relative angular momentum change, 
(\Lt\ — |Lo|)/|io|, with time. Plot (b) shows the evolution of the 
z-componont of the angular momentum, Lz , with time. The data 
is from a Lowered Evans model containing 20,000 particles, 10% 
of which are tree particles. 



case is a Lowered Evans Model with 20000 particles with 
10% Tree particles. 



3.6 Collapse of uniform sphere 

We investigate the performance of the SCFTREE code on 
a system which is not initially in equilibrium. For this pur- 
pose we perform simulations on the collapse of a uniform 
density spherical distribution of particles with random ve- 
locities scaled such that the initial virial ratio of the sys- 
tem |2T/W^Jq = 1/2. Following the thorough investigation 
by Hozumi and Hernquist (1995) of the pure SCF code in 
similar non-equilibrium states, we trace the evolution of 
the virial ratio from its initial value of 1/2. Our purpose 
here is not to perform a detailed investigation of the accu- 
racy of the dynamical evolution but as a qualitative check 
that the combined SCF and Tree codes behave as expected, 
and the virial ratio oscillates about the equilibrium value 
of J2r/Vy|p — 1.0. The results are shown for a typical run 
containing 20000 particles, 10% of which are randomly al- 
located as tree particles. The truncation parameters for the 
SCF part are n = 16 and 1 = 6. 

Figure(H) shows the plot of final density profile of the 
same system after a period of t = 120. The separate profiles 
of the Tree and SCF systems of particles are shown, together 
with the total profile of all of the particles. In the example 
shown the system undergoes homologous collapse (Fillmore 
& Goldreich 1984; Gunn 1977) evolving to a density profile 
of p -^ r"^'^ . Towards the core of the system in this example 
the number of particles is too small to adequately resolve the 
detailed evolution (c/. Hozumi & Hernquist (1995) who used 
hundreds of thousands of particles and obtained adequate 
resolution to resolve the flattening of the core.) 



3.7 Disk + Halo + Bulge 

We have seen from the previous tests that the code is effi- 
cient and accurate even when the tree particles are random 
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Figure 8. Final density profile after i = 120 of the collapsed sub- 
virial (|2T/Vy|o = ^) uniform density sphere. Density profiles are 
individually shown for the SCF and Tree components, shifted 
vertically for comparison. 





Figure 10. This plot shows the absolute separation between two 
non-penetrating systems in circular orbit around each other. Each 
system is an identical Lowered Evans model populated with 2500 
particles. One is assigned to the SCF system, the other to the Tree 
system. This figure demonstrates that the method of expanding 
the SCF potential around its centre of mass every timestep re- 
mains accurate to within 5% over 400 time units. 




Figure 9. The viral ratio, |2T/iy|, evolution of the collapse of a 
subvirial uniform density sphere. The total number of particles is 
20000, with a random 10% being assigned to the Tree system. 



distributed amongst the SCF particles. This is the worst 
case scenario in terms of efficiency. During the force calcu- 
lation on an SCF particle it will form a large list of tree cell 
interactions. Now we move to a case where the tree particles 
form some distinct structure within the SCF system, that of 
a disk embedded within a halo. The study of disk evolution 
in a non-interacting system has been the focus of much re- 
search and requires a good deal more attention than a mere 
subsection to investigate it fully. Here we content ourselves 
with an example of how the varying of A'^haio affects the sta- 
bility of the disk. We construct combined disk-halo-bulge 
models using the method of Dubinski & Kuijken (1995). 
The ratio of masses of the components was disk:bulge:halo 
= 1.00:0.37:12.80, the scale radius of the disc was set to 
Rd = 1.0. The number of disk particles in each case was 
6000, and the number of bulge particles was 2000. In four 
runs the number of halo particles was varied from 13000 to 
10^. Each time the disk particles were assigned to the Tree 
system and the rest to the SCF system. 

Figure hd demonstrates increased stability as num- 



Figure 11. Vertical velocity dispersion, Oz and rms disk height, 
2rms as a function of radial distance at t = 120 . Th e construction 
of the disk— halo— bulge models is described in |3.7| . The three dif- 
ferent runs correspond to the curves (a), (e), and (b) in Figure 

|ii 



ber of halo particles, A^haio is increased. In particular at 
A^haio = 100000 little structure has formed in the disk and 
is unchanged apart from a slight thickening, whereas for 
A^'haio = 13000 the disk is particularly unstable to warping, 
and spiral structure seems to be developing. The thicken- 
ing of the disk shown in Figure [111 These figures confirm 
what we would expect as we change the number of A'^haio in 
a self-consistent simulation. With fewer SCF particles noise 
in the SCF expanded potential will cause fluctuations in the 
halo potential to develop, and thus cause the disk to become 
unstable to warping. 



3.8 Dynamical Friction 

Penetrating encounters between stellar systems provide 
many interesting scenarios which one could fruitfully model 
with our hybrid code. Representing interacting systems sep- 
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Figure 12. This figure plots the Tree particles only after t = 120 to show the effect of changin g th e number of particles in the halo i.e. 
the SCF particles, on the stability of the disk. The construction of the models is described in §BJ. The initial conditions are shown in 
plot (a), with the numbers of particles in the Halo increasing to the right. The top row and bottom row show the x-y plane and the y-z 
plane respectively. 
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Figure 13. The decay in separation between a massive spherical 
primary, treated by the SCF system, and a smaller spheroidal 
satellite, treated by the Tree system. The models u sed for each 



curve are summarised in Table 1. As discussed in §3.8 this fig- 
ure demonstrates the effectiveness of using SCFTREE to model 
dynamical friction. 



arately with the SCF and Tree methods invites the ques- 
tion of how well SCFTREE describes the behaviour of SCF 



and Trie particles interacting with each other in differing 



ways. We address this issue by investigating the ability of 
the code to tackle dynamical friction, specifically by running 
a few simulations of the sinking satellite sort. SCF is used to 
model the primary, and Tree used for the satellite. Although 
the force on each SCF particle due to the Tree particles is 
directly calculated (subject to the tolerance criterion), the 
Tree particles only respond to the SCF particles indirectly 
via the expanded field. One might initially be concerned 
that the lack of direct particle-particle interactions which 
account for the deceleration of the satellite will invalidate 
the use of this code in such cases, however we show that so 
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N^ 



N. 



satellite 



SCF truncation 



(a) Approximate analytical solution valid for large r 

(b) 20000 1000 n = 16 1 = 6 

(c) 20000 1 n = 16 / = 6 

(d) 20000 1000 n = 8 1 = 2 

(e) 20000 1 n = 8 1 = 2 



Rn 



Primary: Hernquist model; B,a = 1 
Satellite: LE model; Ra = 0.1; i?max '■ 
Initial separation on circular orbits =15 



: 3.6: 



= 200: 

M = 



M = 
0.75 



7.39 



Table 1. Summary of the models and model parameters used in 
the tests of dynamical friction 



long as the SCF expansion is truncated after a reasonable 
number of terms dynamical friction is acceptably treated. 

Figure([L3|) demonstrates the dependence of dynamical 
friction on the number of expansion coefficients used in 
the expanded SCF potential. The initial conditions and the 
models used to represent the primary and satellite are sum- 
marised in Table hi Curve (a) in the Figure is a simple ana- 
lytical approximation of the sinking rat e based on the Chan- 



drasekhar dynamical friction formula (Binney & Tremaine 



1987). We approximate /(O) ^constant, and assume instan- 
taneously circular orbits for the satellite. This calculation 
carries with it the assumption of a fixed primary and is 
valid only in the region of large separations, r, low circular 
velocity, no satellite mass loss, and is to be disregarded at 
small r. 

As one might expect when the SCF expansion is trun- 
cated at low orders [curves (d) and (e)] dynamical friction 
is poorly modelled. Increasing the truncation order [curves 
(b) and (c)] with a corresponding increase in the resolution 
of the SCF system causes the Tree system to respond real- 
istically to perturbations in the SCF density field. With the 
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SCF expansion truncated at n = 16 and I — 6 and taking 
a single point mass for the satellite (i.e. one tree particle) 
dynamical friction is modelled well (see curve [c]). We find 
these sinking times are consistent with the results obtained 
by Hernquist and Weinberg (1989) for fully self-consistent 
simulations using a pure tree code and concur with their 
findings [and of White (1983)] that when the response of 
the primary is included in the numerical calculations the 
sinking time increases by a factor of around two. 

Decay of the orbit of a single point mass means in simple 
terms that orbital energy is transferred to dynamical heating 
of the particles in the primary. For an extended satellite 
(i.e. one composed of N satellite bodies), orbital energy can 
also be transferred to the heating of the satellite, thereby 
increasing the decay rate and shortening the decay period. 
This interesting result is demonstrated in curve (b), and 
certainly warrants greater detailed exploration. The result 
has been clearly alluded to in Weinberg (1989) in terms of 
coupling in the analytic linear theory which he describes. 
The extended satellite has many weak resonances which can 
couple at smaller radii, which the point mass does not. 



This implementation of SCFTREE is able to accurately 
model effects of dynamical friction and thus the application 
to systems of sinking satellites is a prime application. Other 
eminently suitable applications will be the evolution of in - 
clined disks in flattened halos (Dubinski & Kuijken 1995), 



the evolution of unstable high-surface-density disks (Dal- 
canton ei al. 1996; Vine & Sigurdsson 1997, in preparation), 
and weak encounters between massive elliptical systems and 
smaller disk systems (Vine 1997, in preparation). 

This preliminary version of the code has much potential 
for future development and expansion, especially in terms of 
performance. The nature of the SCF code makes it suitable 
for parallelisation (Hernquist et al. 1995) together with a 
vectorised version of the Tree code utilising heterogeneous 
systems. With the advent of parallelised Tree implementa- 
tions a doubly parallel SCFTREE is a possibility. 

An interesting development will be to experiment with 
multiple expansions of the SCF systems. With such a tool 
it will be feasible to model interactions between two or 
more large spheroidal systems. For example in the evolu- 



tion of gal axy groups and multiple mergers (Weil & Hern 



Fit ally we note that the sinking satellite loses much of quist 1996). The principle of having more than one expan- 



its mas s during its descent. For the example traced by curve sion centre has been exploited some time ago ( van Albada 



(b) in Figure(13|), 40% of the Tree particles become unbound 



van Gorkom 1977 ) 



from tree system at r ~ 7, and 80% are unbound at r ~ 2. 
The peculiar behaviour of the curves once having reached 
r ~ 1 is due to the practical disruption of the satellite. 



4 DISCUSSION 

In this paper we have presented SCFTREE, a hybrid N-body 
code combining the Hernquist & Ostriker Self-Consistent 
Field code, and the Barnes-Hut hierarchical tree algorithm. 
The S CF code is designed to model svstems with struc- 



Work currently in progress includes a powerful improve- 
ment to the algorithm presented here, that of swapping par- 
ticles between the SCF and Tree system when one system 
becomes more appropriate then the other. For instance the 
disruption of a satellite composed of Tree particles could 
benefit from transfer of escapers to the main body of the 
SCF primary. Conversely, structure formation in A'' > 10* 
cosmological simulations can have increased resolution by 
assigning increasing numbers of Tree particles to structures. 

Another relatively simple development is to incorporate 
hydrodyna mics into the Tree part of the code ( Hernquist & 



ture re jembling to flrst order the Hernquist density proflle, 
P(^) 0*- , _!_ sa ■ The treecode technique is a proven approach 
to effectively model systems that have no a priori and in ef- 
fect, no limit to dynamical or spatial range. The principle 
behind the SCFTREE scheme is to model evolving structures 
that interact with or are within an approximate Hernquist 
potential. The uses of such an approach are wide ranging, es- 
pecially modelling the evolution of structures within a dark 
matter halo, e.g. disks, bars, satellite galaxies, interactions 
and encounters, and cluster galaxies within a cluster halo. 

The prime advantage of this new technique is a signifi- 
cant improvement in performance over using pure treecodes 
which scale as 0{N\ogM). To date only treecode techniques 
have been adequately suited to systems with such large dy- 
namical ranges. However, so much CPU time is expended 
on parts of the system with little dynamical evolution such 
as the halo. By representing the halo with the SCF tech- 
nique, where the CPU time scales as 0{N), we are able to 
reduce the overall CPU time and thus can increase the total 
number of particles. 

As this implementation of the code is based around 
global geometries which behave approximately as B}^'^ in 
projected density, its use would not be appropriate to model 
objects wi th s ubstantially different density profiles. As men- 
tioned in §2.2 other basis expansions for different mass distri- 



Katz 1989). This could be f urther extended to encom pass 
star formation prescriptions (Mihos fc Hernquist 1994). Ul- 



timately it is not unreasonable that the technique could be 
employed in cosmological simulations which are now capable 
of dealing with much more sophistic ated physical processes 
(Katz, Weinberg, & Hernquist 1996) 



butions are possible. However, there are many systems that 
may be modelled by appropriate application of SCFTREE. 
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